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DEMODULATION AND PHASE ESTIMATION OF TWO-DIMENSIONAL 

PATTERNS 

Field of the Invention 

5 The present invention relates to two-dimensional code patterns and more 

particularly to the automatic demodulation of fringe patterns. The present 'invention also 
relates to the automatic estimation of partem and texture orientation of an image. The 
present invention further relates to phase demodulation of a sequence of fringe patterns 
with arbitrary and unknown phase shifts between them. Such patterns typically occur in 
10 optical interferograms, but can also occur in interferograms produced by more general 
waves such as acoustic, synthetic aperture radar (SAR), ultrasound, and x-ray. More 
general image sequences, such as video, can in some cases also be modelled as fringe 
pattern sequences. 

Backgrouad 

Methods exist for demodulating fringe patterns. Examples of such patterns are 
illustrated in Figs. 5 and 14A. However, conventional methods of demodulating closed 
curve fnnge patterns, including circular and elliptical fringe patterns, inev Jt ably produce 
ambigumes. Ambiguities are then typically resolved by including additional constraints 
These additional constraints may include restrictions on the smoothness of curves within 
20 the fringe patterns. 

In addition to ambiguities other defects and errors may also be produced by the 
prehmmary demodulation process. An example of this can be found in the classic Fourier 
(Hilbert) transform method (FTM). When the FTM is applied to a whole image, artefacts 
are produced, the artefacts being related to discontinuities introduced by the half-plane 
25 filter used m Fourier space. The artefacts result in errors in phase estimate, and 
specially m regions where an angle of the fringe pattern is close to perpendicular with 
the half-plane discontinuity line. 

A method adopted to overcome this problem is to choose the half-plane filter in 
Fourier space such that the discontinuity line does not cut through the Slg nal in Fourier 

30 space. 

However, the frequency components derived from circular or elliptical fringe 
patterns are also "crcular" and no discontinuity line could be defined which does not cut 

through the frequency signal. 
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Other methods based on local (small kernel) demodulation may avoid the errors 
of the FTM, but typically fail to disambiguate correctly when there are strong 
perturbations in the underlying fringe pattern. Features such as bifurcations and fringe 
endings represent typically strong perturbations. 
5 Images are often characterised by a number of features including pattern and 

texture. A number of image processing operations depend upon the estimation of the 
orientation of these image features. 

Image processing/analysis applications where the estimation of orientation is 
Q important include: 
S|0 o edge orientation detection; 

j= ° orientation input to steerable filters; 

[: ° pattern analysis; 

13 ° texture analysis, where the local orientation of texture is estimated for 

f 3 characterisation purposes; and 

h 1 5 ° orientation selective filtering, processing, and enhancement of images. 

jL{ 0ne such where the estimation of orientation is particularly important is the 

U demodulation of fringe patterns or equivalently, AM-FM modulated patterns. One area 

where such patterns exist is in fingerprint analysis. Other areas include the analysis of 

naturally formed patterns such as an iris diaphragm, ripples in sand, and fluid convection 
20 in thin layers. The primary reason why it is useful to know the fringe orientation when 

demodulating such a fringe pattern is because it allows for 1 -dimensional demodulators to 

be aligned correctly. 

Another area where the estimation of pattern orientation is used is in image re- 
sampling and enhancement. 

25 Known methods for estimating pattern orientation include gradient-based 

methods. A serious deficiency with such gradient-based methods is that at the ridges and 
valleys of images, the gradient has zero magnitude. Gradient methods are therefore 
unable to provide information about pattern direction in such regions. 

Interferometry refers to experiments involving the interference of two light 

30 waves that have suffered refraction, diffraction or reflection in the interferometer. Such 
experiments typically involve the testing of lenses, mirrors and other optical components. 
The principle parameter of interest in interferometry is the phase difference between 
interfering beams. 
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Present methods of demodulating sequences of phase-related fringe-patterns 
generally rely on a priori knowledge of the phase-shifts between each of the ind.vidual 
patterns. When these relative phase-shifts are available, it is p0SSlble to estimate the 
spatial phase by u Sl ng a generalised phase-shifting algorithm (PSA). In cases where the 
> re lanve phase-shifts are not a priori known, the estimated spatial phase wil, be incorrect 
unless special error-reducing PSAs are used or methods to estimate the actual phase-shift 
between patterns are used. 

Error-reducing PSAs are usefu, where .he deviaSon of ,he ae.ua! phase-shift 
from .he expected phase-shift is sma!, (typicaIlv _ ^ „ , ^ d 

^TT 18 largCT ,Wn - d — « * — *»y randoo,), 

then other methods must be used. 

Methods have been developed «o estimate phase shifts berweei, phas.-rela.ed 

5 c r n of " meuiods - ^ ° n — - °~ ~ — t 

eTT al 5 Pa,Km5 10 ^ «— - *** PW shifts ,o « 

ellipse also require a. ieas, 5 fringe patterns. Mefcods hased on .he s.atis.ica, properUes 

: :i: it :r™ it bcen "~ but - ^ - ^ 

eorre.anon between patterns have been propose, hu, are very sensi,ive „ the „,Ler of 

22Tr m "~ ThMe COI " lafi ° n mah0<15 —~ *« - Ct are 
~ly hiear .„ the frame, which means only ^ « 

patterns typ,caUy reqinre a. leas. 15 frames to operate effeetively 

teJ"^*' ° f of phase-related 

,Hn8e — ' "* ^ *~ * - «* ^ .he case of .hree or 

DiscDosure of the Invention 
It is an object of the present invention to substantially overcome or at least 
amehorate, one or more disadvantages of existing arrangements 

_ According to a fir5t aspect of the invention, there 1S proved a method of 

demodulating a real two-dimensional pattern the methoH , • ■ u 

H ra ' me me thod comprising the steps of: 

pattern & t * n >*™*<«* P*«*™ from said real two-dimensional 

pattern using a two-dimensional spiral phase filter; and >™nsional 
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an. •„ Cre " tlnS " dem ° dUlated ^ ^ -mbining said real two-dimensional pattern 
and said estunated quadrature two-dimensional pattern. 

to -P-t of the invention, there is provided an apparatus for 

demodulatmg a real two-dimensional pattern, the apparatus comprising- . 

5 dan- • T" CStimatin8 3 ^-—^ Pattern from said real two- 

d lmensi onal pattern using a two-dimensional spiral phase filter; and 

means for creating a demodulated image by combining said real two dimen," , 
Pattern and said estimated quadrature two-dimensional pattern ^ *-n S1 ona, 

Accordmg to another aspect of the invention there is provided an a «™- a t * 

~- Po^It^' m0dUla '° r " ^ ^ ^ *- — **« - 

an image sensor for capturing said conjugate of said „,l „ j - 
pattern. J "gate ot said real two-dimensional 

—a, in t:zzz:~zz:iz ~ n ^ - '^ ■ — - 

of . H 311 lma ge, said method composing the steps 

, mage; aPP ' y ' nE * C ° mito ~»' ~ - - taage to provide an energy encoded 

determining a phase component of said energy encoded image, and 
encoded LIT" 8 ^ *~ ~" *~ — - said energy 

ngie or a pattern in an image, said apparatus comprises: 
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means for applying a complex energy operator to the image to provide an energy 
encoded image; 

means for determining a phase component of said energy encoded image; and 
means for calculating said orientation angle from said phase component of said 
5 energy encoded image. 

According to yet another aspect of the invention there is provided a method of 
estimating a spatial phase of fringe pattern images in a sequence of fringe patterns, said 
method comprising the steps of: 

a) removing offsets from each of said fringe pattern images to obtain pure 
10 AMFM patterns; 

b) determining contingent analytic images corresponding to each of said 
AMFM patterns; 

c) determining phase differences from dependent pairs of said contingent 
analytic images; 

15 d) estimating phase shifts between pairs of said contingent analytic images; 

and 

e) estimating a spatial phase of said fringe pattern images. 

According to yet another aspect of the invention there is provided an apparatus for 
estimating relative phase shifts between fringe pattern images in a sequence of phasc- 
20 related fringe patterns, said apparatus comprising: 

means for removing offsets from each of said fringe pattern images to obtain pure 
AMFM patterns; 

means for determining contingent analytic images corresponding to each of said 
AMFM patterns; 

25 means for determining phase differences from dependent pairs of said contingent 

analytic images; and 

means for estimating phase shifts between pairs of said contingent analytic 
images, wherein said phase shifts between pairs of said contingent analytic images are 
said relative phase shifts between fringe pattern images. 

30 Brief Description of the Drawings 

A number of embodiments of the present invention will now be described with 
reference to the drawings and appendix, in which: 

Fig. 1 is an illustration of the local structure of a fringe pattern; 
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Fig. 1; 
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Fig. 2 is an illustration of the Fourier transform of the local fringe pattern of 



Fig. 3 is a flow diagram of a method of phase demodulation; 

Fig. 4 is a flow diagram of an alternative method of phase demodulation; 

Fig. 5 is an example of a fringe pattern image; 

Figs. 6 and 7 are the real and imaginary parts of an intermediate result in 
applying the method of Fig. 4 to the image of Fig. 5; 

Fig. 8 is an alternative estimate of fringe orientation of the fringe pattern image 
of Fig. 5; 

Figs. 9A and 9B are the magnitude and phase components of a demodulated 
fringe pattern when the estimate of fringe orientation of the fringe pattern image as shown 
in Fig. 8 is used; 

Fig. 10 is a schematic block diagram of an optical spiral phase demodulator; 
Figs. 1 1 A and 1 IB are the magnitude and phase components respectively of an 
intermediate result in applying an energy operator to the image of Fig. 5; 

Figs. 11C and 1 ID are the magnitude and phase components respectively of 
another intermediate result in applying the energy operator to the image of Fig. 5; 

Figs. HE and 1 IF are the magnitude and phase components respectively of 
applying the energy operator to the image of Fig. 5 ; 

Fig. 12A is an example of a fringe pattern image containing noise; 
Figs. 12B and 12C are the magnitude and phase components respectively of 
applying the energy operator to the image of Fig. 12 A; 

Figs. 12D and 12E are the magnitude and phase components respectively of 
applying a modified energy operator to the image of Fig. 12A; 

Fig. 1 3 A is another example of a fringe pattern image containing noise; 
Figs. J3B and 13C are the magnitude and phase components respectively of 
applying the energy operator to the image of Fig. 13A, 

Figs. 13D and 13E are the magnitude and phase components respectively of 
applying the modified energy operator to the image of Fig. 13 A; 

Fig. 14A is an example of a frmge pattern image with amplitude and frequency 

modulation; 

Figs. 14B and 14C are the magnitude and phase components respectively of 
applying the energy operator to the image of Fig, 14A; 
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Figs. 14D and 14E are the magnitude and phase components respectively of 
applying the modified energy operator to the image of Fig. 14A; 
Fig. 15A is an example of a texture image; 

Figs. 15B and 15C are the magnitude and phase components respectively of 
applying the energy operator to the image of Fig. 1 5 A; 

Figs. 15D and 15E are the magnitude and phase components respectively of 
applying the modified energy operator to the image of Fig, 15 A; 

Fig. 16 is a flow diagram of a method of estimating the pattern and texture 
orientation of images; 

Fig. 17 shows the Fourier domain coordinates; 

Figs. 18A and 18B show the phase shift parameters S„ for the case where N=5 
on a polar plane; 

Fig. 19 shows a flowchart of an automatic calibration method; 
Fig. 20 shows yet another an example of a fringe pattern; 

Fig. 21 A and 21 B show the magnitude and phase components respectively of an 
analytic image of the example fringe pattern shown in Fig. 20; 

Fig. 22A shows a phase difference between dependent analytic images; 

Fig. 22B shows a histogram of a spread in values of the phase difference shown 
in Fig 22A; 

Fig, 23A shows a modulus of the phase difference between dependent analytic 
images shown in Fig. 22 A; 

Fig. 23B shows a histogram of a spread in values of the modulus of the phase 
difference shown in Fig. 23B; 

Fig. 24 shows a polarity function calculated from the phase difference shown in 
Fig. 22A and the modulus of the phase difference shown in Fig. 23A; 

Figs. 25 A and 25B show estimations of an amplitude modulation and a phase of 
the fringe pattern shown in Fig. 20; 

Fig. 26 shows an estimate of an offset of the fringe pattern shown in Fig. 20; 

Fig. 27 shows a simple weight function calculated from the estimated errors in 
the analytic image shown in Figs. 21 A and 21 B; and 

Fig. 28 shows a schematic block diagram of a general purpose computer upon 
which arrangements described can be practiced. 
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Detailed Description including Best Mode 

Consider a two-dimensional fringe pattern with intensity of the fringe pattern as 

follows: 

f(x.y) = a(x,y) + b(x,y).cos[2n(uo{x-xo}+v 0 {y-yo})+x] ~ (1) 
Position coordinates ( x ,y) may be continuous for an analog pattern or discrete 
for digital patterns. A background level is denoted by a(x,y) while an amplitude 
modulation term is denoted by b(x,y), both assumed to be slowly varying functions- 
Phase term z(x 0 ,y o ) presents the local offset of the overall phase term in the square 
brackets [] of Equation (1). 

Referring to Fig. 1 where a 2-dimensional fringe pattern or AM-FM pattern, 
located around a point (x 0 ,y 0 ) and formed by a number of fringes 10 is shown. A small 
region of interest is defined by £2. Within this small region fi, the spatial frequency 
components u 0 and v 0 are slowly varying so that they are effectively constant around the 
point (x 0 ,y 0 ), making the carrier a linear function of x and y. A normal 12 to the local 
15 fringes 10 is at an angle j3 0 to the x-axis. 

A noise term is omitted from Equation (1), but may be present in real fringe 
patterns due to the occurrence of blurring, non-linearities, quantisation errors, smudging, 
scratches, cuts, dust, etc. 

The intensity pattern f(x,y) is (generally) a rapidly varying function of position 
20 coordinates (x, y) . 

The Fourier transform (FT) of the above fringe pattern region Q is: 

^n(".v)= $lfix,y)ex£-2m{ux+vy)]dxdy 

° (2) 
= Aniu, v) + (B a (u - u 0 , v - v 0 ) exjfor] +B n (u + « 0 ,v+v 0 )exp[-i / fDe- 2rf(, ^* v > i > ) 

wherein i = V— 1 

From Equation (2) it is clear that the Fourier transform F n of the fringe partem/ 
has 3 lobe centres, the first being at the origin and another two at spectral coordinates 
(«o,v 0 ) and (-K 0 ,-v 0 ). The last two lobe centres are mirror images of each other with 
regards to the origin. The effect of windowing the fringe pattern /over a small region Q 
is to blur the otherwise sharp lobes. The central lobe A n> representing the zero frequency 
component, can be removed, and will therefore be disregarded in the following analysis. 



25 
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Methods that may be used for removal of the zero frequency components include using a 
suitably matched high pass filter. The remaining spectral lobes are illustrated in Fig. 2. 

The local spatial frequency components u 0 and v 0 of the fringe pattern/ satisfies 
the following Equations: 



In 

tan(/? 0 )=^ 



(3) 



wherein fa is an angle of the normal to the fringes with the x-axis and ?io is the spacing 
between fringes. 

The conventional Fourier transform method (FTM) of fringe analysis is based on 
the assumption that the angle fa is confined to a small range for the entire fringe pattern, 
typically much less than n radians. However, when considering fringe patterns with 
circular, elliptical carrier patterns or other closed curves, the angje fa covers the entire 
angle range from 0 to 2n radians. The FTM would therefore be deficient in solving fringe 
patterns created with closed curve carriers. The reason for this failure is related lo the 
inability of the FTM to isolate the two (mirrored) lobes of the FT. As noted above, the 
15 frequency components derived from closed curve fringe patterns are also "closed curves". 
Therefore, referring again to Equation (2) s both the second and third terms will have 
"closed curve" spectral components, causing a total lobe overlap. 

The preferred implementation removes directional discontinuities in the Fourier 
processing of the fringe pattem/by smoothly encoding the fringe direction by utilizing a 
filter with odd (180° rotational) symmetry, but (unlike the FTM) a smooth transition 
between adjacent frequencies. Also, the filter used in the preferred implementation does 
not alter the strength of different frequency components, and is therefore a unit magnitude 
or phase only filter. 

A (complex) Fourier representation P( u .v) of the filter preferably has the 
25 following anti-symmetry qualities: 

P(u, v)=exp[iAtw, v)]=-P(-w,-v)=-exp[:(A(w i v)+n)] (4) 
A preferred filter chosen has a phase A(«,v) which is set equal to the polar angle 
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J¥(m,v) = <t> (5) 



where 



u = q cos(^) 
v = #sin(^) 
u 2 +v 2 =q 2 ) 



(6) 



The resulting filter is termed a "spiral phase filter" because the filter phase X 
increases linearly as the polar angle <f> increases. The magnitude of the filter is constant 
(so the filter is a pure phase filter) and the phase has the form of a spiral or helical surface 
with just one turn. An equivalent form of the filter is simply: 



P(u, v) = = 



1 



(7) 



The Fourier transform of the fringe pattern after the central lobe has been 
10 removed is as follows: 

G nM=(B a (u-u 0 ,v-v 0 )ex.^i % ]+B n (u+u 0 ,v+^)ex^-i^^^ (8) 
Applying the spiral phase filter P(u.v) to Equation (8) changes the lobe symmetry to give: 

= 0&n (" ~ "o . v ~ v o ) axpfcflexpp*] -£ n (u + u 0 ,v + v 0 ) exp [-i X ] expgfl^-"*"***** (9) 

Equation (9) can be inverse transformed, as is shown in Equation (10), by 
implicitly assuming that the polar angle ftu.v) « ± A> remains constant over each lobe. 
This is equivalent to assuming that the lobe subtends a small polar angle <f>, which is a 
reasonable assumption for typical fringe patterns/ 

Six, y) m ]" f G(u, v) exp[2m(ux + vy)]dudv 

J_M (10) 
= i exp[i/? 0 ] sin[27r(u 0 {x - x 0 } + v 0 { y - y 0 } ) + x ] 

The desired sine component of the local fringe pattern §(x,y) can then be 
extracted by multiplying g(x, v) with an orientational phase component, exp[-/^b] as 
follows: 

g(x. y).exp{~i/3 0 J =i&in[2a (u 0 {x - x 0 } + v 0 {y - y 0 }) + X J ( n ) 
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Conceptually the process is repeated for all fringe regions with their 
corresponding fringe orientation angles fi. In practice the procedure can be applied to all 
regions simultaneously. It is noted that the original (conceptual) windowing is not 
actually required. 

The aforementioned process depends upon the functions b(x.y), fi(x,y), x ( x ,y), 
u 0 (x,y). v 0 (x,y), and tr&y) having suitably slow variation. However, it is found in 
practice that the preferred implementation can be successfully implemented on any fringe 
pattern/with a suitably smooth variation in orientation J3(x#) and spacing U For fringe 
patterns / with circular fringes and constant or quadratic (chirp) fringe spacing Xo, the 
preferred implementation only breaks down near the very center of the fringe pattern/ 
due to the assumption of small polar angle subtense failing near the centre. 

Equation (1 1) represents an estimation of the imaginary part of the analytic 
function (or im age) related to the original real fringe partem/ This component is also 
called the quadrature component or analytic conjugate. By forming a new (complex) 
output image comprising the original real image/and the estimated imaginary image, the 
phase and magnitude of the resulting complex fringe representation can be obtained as 
follows: 

<f~ *>H§<x,y) exp{-i#>}) m b.exp{i[27n(u 0 {x-x 0 }+ Vo ty-y 0 }y+ % ]y (12) 

Fig. 3 is a flow diagram of a phase demodulation process 100. The fringe pattem 
image/fr*), captured in step 110, typically includes only a real part. The fringe pattern 
image f(x,y) is transformed in step 120 by applying a 2-dimensional FFT. Because the 
fringe pattem image fOcy) is real, the real and imaginary parts of the transform F(u,v) are 
symmetric and antisymmetric respectively. 

The DC component (lobe) of the transform F(u. v ) is removed in step 130. In 
step 40, the spiral phase filter P(u,v) is applied to the output from step 130. In step 150, 
the inverse. Fourier transform is calculated. This may be performed by again applying a 
2-dimensional FFT. 

In step 160, an orientational phase filter is applied to the resulting image of step 
150. In step 170 the real part is discarded, leaving only the imaginary part, which is then 
added to the input fringe pattem image f(x,y) in step 180 to produce an output image in 
step 190. 

Fig. 4 is a fl ow diagram of ^ alternative i mp i ernentatjon of a phage 
demodulation process 200, and will be explained with reference to an example fringe 
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pattern. An example of such a two-dimensional fringe pattern is illustrated in Fig. 5, 
wherein only the real part of the pattern is known. 

The fringe pattern image f(x.y) is captured in step 210. In the phase 
demodulation process 200, steps 120 to 150 of the phase demodulation process 100 
shown in Fig. 3 are replaced by a single convolution step 220. The fringe pattern image 
f(x.y) is preferably convolved with spiral phase function p(x,y), wherein p(x,y) is denned 
as follows: 

P(*.y) = £ £/»(«; v)exp[2**(ux + vy)]dudv 

£° r> (u + i'v) r „ „ (13) 

P ( x ,y) = -L. — ZL = exp(ie) 



wherein 



x = r cos 0") 

y = r sin 0J (15) 

The result after such a convolution on the fringe pattern image f(x.y) of Fig. 5 is 
illustrated in Figs. 6 and 7, showing the real and imaginary parts otf(x,y)*p( x ,y). 

Referring again to Fig. 4, in step 230 the orientational phase filter, exp[-ip] is 
applied. However, an estimation of the fringe orientation p is needed. 

One method of estimating the fringe orientation p is to compute the 2-D energy 
operator which gives uniform estimates over the fringe pattern image f(x,y). In contrast, 
typical gradient based methods fail to give reliable orientation estimates at the peaks and 
valleys of the fringes. 

An alternative method of estimating the fringe orientation p is to use an uniform 
orientation estimator, but to utilise local environmental constraints (such as topology of 
smoothly varying fringes) to resolve or unravel the orientational ambiguity. Fig. 8 
illustrates a smooth estimate of the fringe orientation exp(</7) derived from exp(2^) . The 
orientation phase estimate exp(,/?) has only one discontinuity line at positions where the 
phase changes from 360' to 0°. A number of different methods are suitable for resolving 
the orientational ambiguity. The magnitude and phase components of the result of step 
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230 using the estimation of the fringe orientation |J of Fig. 8 is illustrated in Figs 9 A and 
9B. 

It is noted that an estimate of the orientation can be obtained from the output of 
step 150 in Fig. 3 or alternatively from the output of step 220 in Fig._.4. However, 
5 additional unravelling will be required to obtain a useable orientation estimate, 

A preferred implementation of estimating the fringe orientation (J -is described in 
detail below. 

In step 240 the real part is discarded, leaving only the imaginary part, which is 
then added to the input fringe pattern image/f*,^, in step 250 to produce an output image 
10 in step 260. 

The method of Fig, 3 can be practiced using an optical spiral phase demodulator 
300, such as that shown in Fig. 10. The optical spiral phase demodulator 300 is 
illuminated by coherent light 301. A first spatial light modulator (SLM) 302 modulates 
the coherent light 301 to produce the input image, which is typically a fringe pattern 
ffcy). The fringe pattern/^; is Fourier transformed by lens 303. A second SLM 304 or 
a spiral phase plate performs phase modulation in the image projected on it. The phase 
modulated image is then Fourier transformed by lens 305 and is projected onto a final 
SLM 306 or spiral phase plate. The final spiral phase plate 306 retards the coherent 
optical beam in a manner such that the plate introduces a spiral phase multiplication. An 
image sensor 307 is in close proximity to the final SLM 306 to detect an output image. 
The output image will be in quadrature (or analytic conjugate) to the input image. In 
other words, if a cosine pattern is input, a sine pattern is output. 

In the special case where the input image is a circularly symmetric fringe pattern, 
then the SLM 306 has the simple form of a spiral phase with opposite sign to SLM 304. 
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Returning to the description of the estimation of the fringe orientation p, it is 
noted that in recent years, an efficient method of signal demodulation using a concept 
known as an "energy operator" has become popular. The operator can be used for 
demodulation of AM and FM signals, including combined AM-FM signals, using a result 
from the analysis of simple harmonic motion (SHM). The energy operator ¥ is defined 



as: 



533242US.doc 



10 



- 14- 

A particularly usefbl property of the energy operator »P is that it gives uniform 
estimates of modulation parameters anywhere on the signal / Essentially this is due to 
the fact that, at local peaks (and troughs) of the signal /, the energy is concentrated in the 
potential energy component, represented by the second term in Equation (16), whereas at 
the zero crossings, the energy is in the kinetic component represented by the first term in 
Equation (16). In other words, the rapidly varying sine and cosine components associated 
with AM-FM signals are exactly cancelled out by the energy operator leaving an 
underlying constant representing total energy. 

The operation of the energy operator ¥ has also been extended to N-dimensions. 
This is described in the publication "A multi-dimensional energy operator for image 
processing", SPIE Conference on Visual Communications and Image Processing, Boston, 
MA, [1992], 177-186, by P. Maragos, A.C. Bovik, and J. F. Quartieri. 

Utilising the energy operator ¥ on images, (i.e. data sets formed in 2- 
dimensions.) separate x and y components may be calculated and added together to obtain 
1 5 an overall 2-dimensional energy operator. 

However, the energy operator ¥ defined in this simple way has certain 
deficjencies. The main deficiency being that the orientation associated with the energy 
components is strictly limited to the positive quadrant of the x-y plane, which means that 
visually very different orientations are indistinguishable using two positive energy 
20 components. 

Referring again to Fig. 1 and Equation (1), the intensity of the fringe pattern can 
be represented as follows: 

gix.y) = b(x,yj.cos[2n(u 0 {x-xo} +v 0 {y-yo} )+z] ( 1 7) 

wherein any slowly varying (non-multiplicative) background offsets a(x,y) are removed 
with an appropriate high-pass filter. The local spatial frequency components u» and v 0 
satisfy Equations (3). 

The x and y component energy operators ¥ x and % respectively may be 
performed on the fringe pattern defined in Equation ( 1 7) as follows: 



25 



& -s-^r-ww (18) 
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Substituting Equations (18) and (19) into Equations (3), the local frequency and 
orientation of the local fringe pattern may alternatively be expressed as follows: 

(2xbf (20) 



¥1 



to < 21 > 

Equation (21) indicates a loss in sign information, which corresponds to an 
inability to distinguish between fringes at four different orientations 
Po--Po>x + P<»x-Po- This is very undesirable because the human eye can clearly 
separate orientation 0 O from and - A from n + fi 0 . In general, however, it is 

not possible to distinguish (on a local level at least) between any orientation and the 
it rotated orientation. 

An alternative approach to extending the energy operator V to two dimensions is 
# by defining a two dimensional complex operator T c which encodes the energy in the 
magnitude and the orientation in the phase: 

" vA s ) + v y {g} = {z>*Tte+ivl) (22) 

As mentioned above, the simple x and y components alone fail to represent the 
orientation correctly. A preferable approach would be to have a 2D energy operator that 

has the following output: 

V c {g(x)} = (2nbJ (u 0 + | Vo ) 2 - {2yra 0 bf exp(2,/? 0 ) (23) 

The complex energy operator * c enables the isolation of the fringe spacing Xo, as 
seen in Fig. I. and orientation angle p 0 respectively. It is noted that the term 2 in the 
exponential factor of Equation (23) again causes the orientation angle 3o to be defined 
modulo n. Such an energy operator is developed by defining a 2- dimensional complex 
differential operator D as follows: 

XSS dx By (24) 
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The complex energy operator T c can then be defined as follows: 

( D k}y-g^{g} = (2^ o bfcxp(20 o )=^ e { g } (25) 
With the Fourier transform defined as follows: 

G (u, v)= J jg( x , y)exp[- 2m(ux + vyypxdy (26) 

the operator D can be implemented in the Fourier space , providing the following 

relationship: 

CO «o 

(2*2X" + 'v)G( M ,v)= //Z)fe(x^)}exp[-2/n( l « + v V )}i e ^ (2 7) 

— <3J — oft 

wherein 

2;r(u + iv) = 2^. exp(i>) 
and the polar coordinate transforms are: 



(28) 



u—q cos(#>) 
v=^sin(^?) 
M 2 + v 2 =^ 



(29) 



The corresponding Fourier operation of Equation (27) provides for the 
multiplication of the Fourier transform of the fringe pattern G(u,v) with a spiral phase 
filter, tfexpM The spiral phase filter has the characteristic of linearly (conical) 

increasmg in magnitude with frequency. The Fourier transformation pairs may be written 

as: 

(2m)qexp(i<p)G(u,v)<z>D{g(x,y)} (31) 
(2m) 3 9 2 exp(2 l >)C7(« > v)^> J D^ g . ( ^ > , ) } (32) 

The derivatives in Equation (25) can therefore be implemented as 2-dimensional 
Fourier filters with spiral phase and conical magnitude. 

By applying the differential operators defined in Equations (30) to (32) as filters 
through simple multiphcation of the image of Fig. 5, the results as shown in Figs 1 1 A to 
HF is obtamed. Figs. 11A and 11B show the magnitude ^ phase componems 
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respectively of the {D{g)f term of Equation (25), calculated by using the relationship in 
Equation (31). Similarly, Figs. 11C and UD show the magnitude and phase components 
respectively of the gD 2 {g} term of Equation (25), calculated by using the relationship in 
Equation (32). The magnitude and phase components respectively of the complex energy 
operator ¥ c is shown in Figs. 1 IE and I1F. Referring to Equation (23), the phase 
component of the complex energy operator T c has double the orientation angle j3 0 of the 
fringe pattern / 

All the above calculations and derivations may alternatively be implemented as 
convolution operators without Fourier transformation. 

However, when the energy operator ¥ c is applied to a fringe pattern /containing 
noise, such as that shown in Fig. 12A, the resulting magnitude and phase components of 
the complex energy operator *f 0 , shown in Figs. 12B and 12C, are less satisfactory. The 
primary reason for the degradation is the linear (conical) spectral enhancement of the 
differential operator D. The differential operator D effectively enhances high frequency 
noise, as it essentially operates as a high-pass filter. In particular, in the areas of low 
spatial frequency <r 0 , such as the area 20 on Fig. 12A, the contribution by the noise signal 
due to the differentiation dominates the determination of the orientation of the fringe 
pattern. Referring to Fig. 12C, the phase component, from which the orientation can be 
obtained by using the relationship from Equation (23), provides near random noise in the 
corresponding area 22. 

One possible correction is to apply a low pass filter function to the input image 
or the result. However, the energy operator ¥ is not a linear (shift invariant) operator. It 
is dependent on products of the input function g and its derivatives dg/dx and tfg/dx 2 , as 
can be seen from Equation (16). Low pass filtering results in the image becoming 
blurred. 

To overcome the deficiencies of the complex energy operator T c , while 
maintaining the ease of calculation of the orientation angle fy. a pure phase operator D M , 
thus having uniform (spectral) magnitude, is created by setting q in Equations (31) and 
(35) to unity, to define the following Fourier transform pair: 

ex.p(i<p)G <=> D M {g(x, y)} (33) 
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The pure phase operator D M is a spiral phase (or vortex) operator which results 
in the removal of noise without significant blurring of the signal. By applying the 
modified energy operator defined as: 

**{g}~<Pu{g)Y-gDl{s} - (34) 

to the image in Fig. 12A, the magnitude and phase components respectively of the 
modified energy operator T M shown in Figs. 12D and 12E, are obtained. The comparison 
of the phase components of the complex energy operator ¥ c shown in Fig. 12C and 
modified energy operator V M shown in Fig. 12E, and in particular to the corresponding 
areas 22 and 24 respectively, indicate that the modified energy operator *P M provides a 
better result. 

It is noted that, similar to the energy operator the modified energy operator 
V M is also non-linear. A consequence of this is that it is generally not possible to obtain 
the performance of the modified energy operator »F M by merely applying a low pass filter 
to the input function, nor by filtering the output function, of the energy operator V. 

The modified energy operator W M has a property of scale invariance, which it 
inherits from the scale invariant £> M operator. Scale invariance is another aspect of the 
operator's spectral neutrality and means that the operator essentially gives outputs 
independent of feature size. 

To further illustrate the advantages of the modified energy operator T M , the 
complex energy operator *F C and the modified energy operator V M are respectively 
applied to an image with uniform noise in the range ±64 greylevels shown in Fig. 13 A. 
The magnitude and phase components respectively of the complex energy operator ¥ c are 
shown in Figs. 13B and 13C, while the magnitude and phase components respectively of 
the modified energy operator T M are shown in Figs. 13D and 13E. Because of the high 
level of noise in the input image, the phase component of the complex energy operator V c 
is completely weighed down by the contribution of the noise. The modified energy 
operator ^ M gives a much clearer phase without any obvious blurring that would be 
expected of noise reduction by conventional low pass filtering of the complex energy 
operator 

Another comparative example wherein the performance of the modified energy 
operator * M can be seen is shown in Figs. 14A to 14E. Fig. 14A shows an amplitude and 
frequency modulated image. The magnitude and phase components respectively of the 
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complex energy operator ¥ c are shown in Figs. 14B and 14C, while the magnitude and 
phase components respectively of the modified energy operator T M are shown in Figs. 
14D and 14E. The complex energy operator ^ c again has difficulty in regions of low 
spatial frequency Xq, such as area 30, causing the corresponding area 32 in Fig. 14C to 
have high levels of uncertainty. On the other hand, the corresponding area 34 when the 
modified energy operator is applied, provides a much clearer result. 

It is further noted that since the underlying fringe pattern of the images in Figs. 
12 A, 13A and 14A are the same with different processes applied thereto, the processes 
not affecting the orientation of the fringes of the underlying pattern, the results of the 
modified energy operator »F M shown in Figs. 12E, 13E and 14E are substantially similar. 

The modified energy operator may be implemented entirely in the Fourier 
space domain by using complex convolution kernels that correspond to the pure phase 
spirals. With the kernel k and its Fourier transform K, the Fourier transform pair can be 
defined as: 

«fr..).«^rt-ji±s: < .i^l.ia^. 4fc ^ (35) 

wherein the spatial polar coordinates satisfy the following: 
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20 



x 2 +y> =r 2 ' 
x = r cos(#) 
y = rsin(0) 



(36) 



The convolution kernel k reduces by the inverse square of the radius, and also 
has a spiral phase structure, but rotated by 90° with respect of the x-y coordinates. 
According to the above, since it is now known those operations required to be performed 
to extract local orientation information from a fringe pattern, it is not necessary to 
perform such calculations in the Fourier domain. Rather, the desired functions may be 
recalculated in the real domain thereby facilitating ease of processing. In many cases the 
desired functions may be approximated by suitably truncated kernels. The inverse square 
25 drop-off allowing relatively small truncation errors even for small kernels 

Referring again to Fig. 1, by considering a sufficiently small region Q of images 
composed of line or edge features, the local structure fits the model of Equation (1). This 
is true for many image features, the most obvious exceptions being textures that contain 
multiple AM-.FM modulated signals. Sometimes an image may approximate the mode] of 
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Equation (17) except for a low or zero frequency background offset in intensity. In such 
cases this background level may be removed by high pass filtering (or related methods) 
leaving an image that does satisfy Equation (17) which can then be processed using the 
complex energy operator formalism. 

Where Equation (17) is representative, the modified energy operator *F M , as 
defined in Equation (34), may be used for the calculation of the orientation -angle fa of all 
such images. 

An example of such an image and the operation of the modified energy operator 
T M are shown in Figs. 15A to 15E. Fig. 15A shows a texture image. The texture image 
is characterised by distinct lines 40 and areas of low variation 41. The magnitude and 
phase components respectively of the complex energy operator ^ c are shown in Figs. 1 5B 
and 15C, while the magnitude and phase components respectively of the modified energy 
operator ^ M are shown in Figs. 15D and 1SE. The phase component of the modified 
energy operator ¥ M provides a result that is more intuitively representative of the 
15 orientation of the texture image of Fig. 15A. 

A flow diagram of a method 400 of estimating the pattern and texture orientation 
of images is shown in Fig. 16. An input image is inserted in step 405. The input image is 
typically digital format with a predetermined number of columns and rows. Alternatively 
the image may be scanned to convert it into the digital format. 

Although the preceding analysis has been in terms of continuous functions and 
. continuous operations, the analysis is also directly applicable to discretely sampled 
functions (images) and discrete operations. 

A Fourier transform of the image is calculated in step 410, typically using the 
Fast Fourier Transform (FFT). The Fourier transformed image is multiplied with a term 
25 exp(i<p) in step 420. 

An output of step 420 is again multiplied with the term exp(i<p) in step 430 
before it is inverse Fourier transformed, typically using an Inverse Fast Fourier Transform 
(IFFT) in step 440. The result is then multiplied with the input signal in step 470. 

The output from step 420 is also inverse Fourier transformed in step 450 before it 
is squared in step 460. The output from step 470 is subtracted from the output of step 460 
in step 480. From . this complex signal, the angle of the signal is calculated in step 490 
before it is halved in step 494 to provide an orientation angle fa as an output in step 496. 
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As noted above, the orientation estimate is modulo n but can be unwrapped using 
continuity constraints to obtain an estimate that is modulo 2n. 

Interferometry is the use of interference phenomena for measurement purposes. 
This allows for the measurement of very small angles or small displacements of objects 
relative to one another. An interferometer is a device to make such measurements, which 
typically uses a beam of light, coming from a single source, such as a star, a laser or a 
lamp, and two or more flat mirrors to split off different light beams. These beams are 
then combined so as to 'interfere' with each other. The beams interfere with each other 
by constructively adding together and cancelling each other out, thereby forming 
alternating bands of light and dark, called fringe patterns. 

What makes the interferometer such a precise measuring instrument is that these 
fringes are only one light-wavelength apart, which for visible light corresponds to about 
590 nm. With the fringe patterns having a cyclic nature with phases that shifted, the 
important measurement parameter is the relative phase-shifts between these phase-related 
fringe patterns. 

Consider an intensity f n {x,y) of the n* fringe pattern in a sequence of phase- 
related fringe patterns with the following mathematical form: 

/* (*. y) = a (x, y) + b(x, y)cos[z(x, y)+S n ] (3 7) 

wherein a{x,y) is again denotes a background or offset, b(x,y) the amplitude 
modulation, and %(x. y) the phase. Each fringe pattern has a phase shift parameter S„ . An 
example of such a fringe pattern is shown in Fig. 20. 

An improved method of estimating relative phase-shifts between phase-related 
fringe patterns/, is proposed. The method is based on an isotropic quadrature transform. 
The isotropic quadrature transform is used to estimate the spatial phase [xix,y)+6 n ] in 
one fringe pattern/,, followed by comparing the spatial phase [x(x.y)+d n ] between fringe 
patterns /„ to estimate the relative phase-shift(s). With 3 or more frames, the method is 
independent of background variations a(x,y) . 

The exact value of phase z(x,y) can be calculated from three or more fringe 
patterns/, if the values of the phase shift parameter S„ is known for each fringe pattern/,. 
For more than three fringe patterns /,(*. y) , the phase x(x,y) is over-determined and is 
usually estimated in a maximum likelihood process. Other estimation processes may be 
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more appropriate in cases where additional systematic error or distortion effects apply to 
Equation (37). 

Consider a least squares estimate for the phase %(pc,y) from N fringe patterns: 

'p 5>. 5>. Y a ) fX/.) 

Z c « Z c » Z^ c « + X = 2>*/- (38) 

.Z j » E^ c n Z j « X- bco *z) lZ*«/«> 

where c„ =cos<5,,, and = sin<5' ) ,, and the summations are over N fringe patterns ie. 

Z = Zl, • squares Equation (38) can be inverted, providing c„ and ,„ are such 

that a singular matrix does not occur, as follows. 



+ 6cos^ 
1^-6 cos X ) 



'Z* 5>. Z-. 
IX Z c * Z^» 

,Z 5 » Z 5 « c .. z 



V 1 


f Z/- ) 




Z c «/» 







(39) 



Equation (39) allows for the offset a(x,y), the amplitude modulation b(x,y), 
and the phase xix.y) to be calculated. However, before Equation (39) can be used, all the 
phase shift parameters S„ must be known. 

The phase shift parameter S„ of each fringe pattern is estimated by using a 2-D 
quadrature technique based on a spiral phase Fourier operator. Estimation of phase xix.y) 
in the above manner has many advantages related to the extraordinary demodulation 
properties of the spiral phase Fourier operator. 

A first step in determining the phase shift parameters S„ is to remove the 
background (offset) a(x,y) from each of the fringe patterns f m (x t y). In the preferred 
implementation, this is done by forming inter-frame differences between pairs of fringe 
patterns f n (x,y) as follows: 

8**,= f„ -/„ = -g mn = b{cos[z + cos[% + <?„]} (40) 

The inter-frame difference function g nm is a pure AMFM (amplitude modulation 
and frequency modulation) function without offset and therefore can be demodulated if 
the inter-frame difference function gnm , which is also a fringe pattern, is locally simple. 
For the fringe pattern g nm to be locally simple, it has to have a unique orientation 
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/3(x,y) and spacing A(x,y) at each location (x,y), except, perhaps, at a finite number 
of locations where the orientation 0(x,y) and/or spacing A( x ,y) may be uncertain. 
Advantageously, local simplicity is a property which distinguishes fringe patterns from 
other, more general, patterns. Fig. 1 shows the local structure of a typical fringe pattern 
for the small region £2. 

A spiral phase Fourier operator V{}, hereinafter referred to as a vortex operator, 
is defined as follows: 

v {/ (*. y)} - F" 1 {exp [,*]F {/ (x, y)}} 

*{f(x>y)} = F(u,v)= j j f{x,y)eKp\-2jzi(ux+vy)~]dxdy 

u = qcos$\ 
v = q sin <p J 

The Fourier operator F{ } transforms a general function f(x,y) from the spatial 
domain ( x , y ) to the spatial frequency domain (u,v). The Inverse Fourier operator 
F"' { } transforms a spatial frequency domain signal back to the spatial domain (x,y). 
In practice, for discretely sampled images, the Fourier and the Inverse Fourier transforms 
are efficiently performed by a Fast Fourier transform algorithm and an Inverse Fast 
Fourier transform algorithm. 

The spiral phase exp [if] is wholly a function of the angular component of the 
spatial frequency polar coordinates ( ff .*) shown in Fig. 17. The vortex operator V{} is 
applied to the inter-frame difference functions g nm to obtain: 

v{ gmn } £ t exp[i^>[ s inOir + <?„)- sinGr + 6 m )] (42) 

From Equation (42) it can be seen that the vortex operator V{} is an isotropic 
quadrature operator in that it converts cosine fringes at any orientation 0( x ,y) into sine 
fringes at the same orientation 0{x,y). The approximation in Equation (42) is-valid 
except in regions where the amplitude modulation b(x,y) , spacing A(x,y) or orientation 
P{x,y) vary too rapidly. Areas where the radius of curvature of the fringe pattern g 
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smaller than the fringe spacing A(x,y) also give a lower accuracy approximation. 
However, for most typical fringe patterns g nm , the approximation is good and any errors 
are isotropically distributed. In contrast, conventional methods, such as the half-plane 
Hilbert transform, cause highly directional errors, and are not suitable for'high precision 
demodulation of general fringe patterns. 

Next, the local orientation p (x,y) is compensating for. To do so, the orientation 
P (x,y) must be estimated. There are a number of ways to do this. One method relies on 
the gradient of the fringe pattern g mn . The ratio of the x and y components gives the 
tangent angle /3(x,y) from the local Taylor's series expansion of the phase X {x,y) about 
the point (x 0 ,y 0 ) : 



X{x, y)= Zoo + (*- *o ) + 2*v 0 (y - y 0 )+ o(r 7 ) 



dx 
dz(x,y) 



= 2nu 0 

jr=jr 0 



(43) 



y*>y 0 



s2m/ 0 



The gradient components are: 

s b{- sinOir + 5 n )+sin( z + Sm ))& 
^M-s m (* + 0+sm(, + 0)|^ ^ 
The ratio of the gradient components gives the tangent angle: 

By I dx "« 0 - tan ^- (45) 

Again the approximation of Equation (45) is valid in regions with smoothly 
varying parameters as defined for the vortex operator V{}. Other methods of orientation 
estimation may be used. For example, methods based on the multi-dimensional energy 
operator or the vortex operator itself may be used. Generic orientation estimation is 
assumed at this point, and the orientation fi for any fringe pattern /„ or inter-frame 
difference function g„ m may be calculated. Alternatively the orientation 0 may be 
estimated from combinations of the previous methods in a statistically advantage manner. 
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It will be assumed hereinafter that there is just one orientation estimate, which shall be 
called orientation estimate fi.. The orientation estimate fi, typically contains an 
ambiguity in the sign: 



exp[//J] = cosjff +/sinyff = + u o + iv ° 



(46) 



25 



Different orientation estimates fi, typically have the sign ambiguities in 
different places. It is worth noting that the overall sign ambiguity arises from the general 
difficulty in distinguishing a fringe at orientation from a fringe at orientation 0 + x . 

Equation (40) is rewritten as: 

g nm =26sin[^ + (^-f^)/2] s in[(^-c?J/2] (47) 
and Equation (42) as: 

exp[- ifi. Ms M1 } s 2ib exp[,(/? - p t )] COS [^ + (S m + 5„ )/2]sin[(<5 w - S n )/ 2 ] (48) 

The complex exponential in Equation (48) has possible values {-1;1}, depending 
on the orientation estimate fi, , and shall be called a polarity factor h(x,y) : 

*(*»jO = «p[>(/»-A)] (49) 

A contingent analytic image is defined as: 

S„ m - i(g nm - exp[- ip c \j{g m }) (50) 

By substituting Equations (47), (48) and (49) into Equation (50), the contingent 
analytic image g nm is rewritten as: 

=2bsin[(S m -S n )/2\h COS \ x + {S m + 5„)/ 2]+ i sm[x + {S n +S n )/l]) (51) 
This definition of an analytic image is contingent upon the orientation function 
A used in the definition shown in Equation (50). It is noted that the contingent analytic 
image g nm can be achieved by placing the inter-frame difference function g nm in the 

imaginary part of a complex image, and placing the imaginary output of the vortex 
operator V{} in the real part. 

A phase a nm of the contingent analytic image g nm is calculated as follows: 

*nn, = Arg{ig nm -iexp[-i0 e ]V {*„„,}}= h\x + {Sm+Sn)/2\ ( 52 ) 
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The of phase ^ is possjMe ^ 1(mg _ ^ _^ )/2>o whi<;h 

only occurs if 2 ^ whe „ y js _ ^ ^ ^ ^ ^ ^ 

«_ » such . cas. „ due „ phase shife ^ ^ ^ ^ ^ 
inter-frame difference function g m being identically zero. 

The phase difference between dependent pairs of contingent analytic images 
such as g nm and ^ , is calculated as follows: 

- a "> =^ + ^«. +6MA-h\x + { 5k + 5J/2J=*„ - 5k )/2] (53) 
This allows a difference i„ phase between any two frmge patterns /„ and/, to be 
found, namely {Sn -S kX with a sign ambiguity as a result of the polarity factor h. It is noted 
*at al, phase calculations, including phase differences/are calculated modulo 2* 
Var,ou s methods may be used to extract both the polarity factor * and the phase 
differences (<?„-<%). In a simplest implementation, from Equation (53), the absolute values 

ot the phase difference between two fringe natter™ f an * r « 

ge P aite rns/„ and f k are averaged over the full 

fringe pattern area S: 



jj s dxdy (54) 

A PrcfeiTCd i-P^entation uses a weighted integral with weighting „( x ,yJ over 

the frmge pattern area A The weightmg ^ used should be a measure of the estimated 

accuracy of the local estimate of the phase cu Thi« ic i . ^ 

pndsc ow This is related to the AMFM 

characteristics of the fringe patterns Ux , y) in the sequence of phase . related 
Patterns. Thus, the absolute values of the phase difference (W0 between two fringe 
patterns/, and/,, weighted averaged over the full fringe pattern area S is: 

Sly(x,y)dxdy C55) 

Th.s allows for a very low (or zero) weighting „( Xl y) to be used where the phase 
a varies quickly, as is the case for d i SContinuities> Qr ^ ^ ^ 
y) is low. Equation (53) can be rewritten as: 
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By arbitrarily taking the phase difference for these particular values of n 

and k as positive, from Equation (56) follows an estimated polarity factor h c as: 

h c (x,y)= "™~ a "* { (57) 

This estimated polarity factor function h e can be used in all further computations, 
5 because a global sign convention has effectively been set for the phase difference (S^-Sk) 
in the fringe pattern sequence. 

Next, the phase difference S n —S k may be calculated for all possible pairs. The 
minimum requirement is to calculate the phase difference S n — 6 k for all sequential pairs, 
which amounts to N values for N fringe patterns f„. By also calculating the phase 
10 difference S„ — S k for non-sequential pairs, additional statistical robustness in the 

estimates is achieved. The number of combinations of two-phase differences S rt ~S k 
from N phases are: 



2(N-2)\ 



(58) 



In the particularly interesting case where N=3 S the number of combination is also 
15 3, which equals the number of sequential differences too. Figs. 18A and 18B show the 
phase shift parameters S„ for the case where N— 5 on a polar plane. The points 1-5 are 
thus on a unity circle at points corresponding to &xp(i5 n ) The connecting line between 
points n and m corresponds to [exp(z<S„)-exp(z <?„)]. The length of each connecting line is 
proportional to sin£(<5 m — 5 n )/2^, which also appeared in the inter- frame difference 
20 function g nm of Equation (47). The length of each connecting line is a direct indication 
of the reliability of the estimate of the phase difference S n —S k , as the length is the factor 
appearing in the magnitude of the inter- frame difference function g nm . Fig. 18A shows 
the connecting lines for all sequential pairs, whereas Fig. 18B shows all possible 
connecting lines, and thus inter- frame differences S n -S kJ which adds another 5 phase 
25 differences S n —5 k - With the phase differences y mn defined as y mn = 5 m - 5 n , an 
estimation y mn of each of the phase differences y mn may be calculated using a weighted 
mean, and compensated by the estimated polarity h e as follows: 
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jj^dxdy 



(59) 

mod2jr 



In situations where the fringe patterns f„ contain noise and/or other distortions, 
the best phase differences estimates may come from a two dimensional fit of all the 
points corresponding to exp(iA), as shown in Figs. 18A and 18B, to a unit radius circle. 
5 Without loss of generality, one of the phase shift parameters S n may be set to 

zero. For convenience the first phase shift parameter S i = 0 is set. This is consistent 
with the idea that the absolute value of the phase %(x,y) is undefined and only phase 
differences =S m - S r are really meaningful. The remaining (relative) phase shift 
parameters S n are calculated as follows: 

S } =0 

These values of the phase shift parameters S n are substituted back into the 
general PSA defined by Equation (39). This allows for all the fringe pattern parameters, 
namely the offset a(x,y), the amplitude modulation b(x,y), and a phase x\x,y) to be 
calculated. It is noted that the calculated phase X \x,y) then includes a common phase 
1 5 shift, namely the phase shift S\ . 

However, there may be some residual calibration error due to the original vortex 
operator failing in certain areas. These errors are rectified by re-evaluating the weight 
function w(x,y) based on the present estimates of the offset a(x t y), the amplitude 
modulation b(x,y) , and phase Z '(x,y) , as well as their partial derivatives. This allows 
for the assessment of how well the original fringe patterns, in the form of the inter-frame 
difference functions gnm , fulfilled the smoothness requirements of the vortex operator 
V{}. Areas where the phase a„ m varies too quickly, or areas where the modulation 
b{x, y) is too low, the weighting function *>(x,y) is reduced or set to zero. Once a new 
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weight function w(x,y) is formed, the phase calibration process, from Equation (55) 
onwards, is repeated, resulting in better estimates for the phase shift parameters S n . 

The accuracy of the estimates for the phase shift parameters S n may be 
expressed as the variance of phase difference, as follows: 

\2 



- , Sis (K" ~ a 'nk\-a nk ) H- 2 (x, y) dxdy 
Sj s ™ 2 (x,y)dxdy 



wherein 



JJ , Kh-^K (x 9 y)d xdy 
jj s *> 2 (x,y)dxdy 



= f f x". / (62) 



Care is needed in evaluating Equations (61) and (62), as the phase xX x *y) is 
periodic and false discontinuities can disrupt the calculations if 5^ near - 7t or n . 
io Other methods for estimating these parameters a„ k that specifically incorporate 

the cyclic nature may be used instead of Equation (62), for example: 

&njc =arg{J^exp[z|tt m „ -a mt |]w 2 (x,y)dxdy} (63) 

Fig. 19 shows a preferred embodiment of a automatic phase-shift, calibration 
method 600. A first step 610 receives a sequence of fringe patterns f n . Each of the 

15 fringe patterns f n is real valued. An example of one of the fringe patterns f n from the 
sequence is shown in Fig. 20. Step 620 removes the fringe pattern offset a(x,y) from 
each of the fringe patterns f n . This is done by forming inter-frame differences g nm as 
defined in Equation (40). The inter-frame differences g nm are also real valued. 

Step 630 forms analytic images g nm corresponding to the inter-frame differences 

20 g nm as defined in Equation (51). This is preferably achieved by placing the inter-frame 
difference function g„ m in the imaginary part of a complex image, and placing the 
imaginary output of the vortex operator V{g nm } in the real part. Fig. 21 A and 2 IB show 
the magnitude and phase components respectively of an analytic image g nm of the 
example fringe pattern f„ shown in Fig. 20. To enable the values of the phase components 
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of the analytic image to be displayed, the values have been adapted to a range 

[0;255] with 0 representing -k and 255 representing — n . 

128 

Step 140 extracts the phase difference a^-cc^ (modulo 2k) from two 
dependent analytic .mages g hm and g niA . Fig. 22A shows a phase difference a -a 
between dependent ana.yt,e image, g „ md f- of the example fringe pattern /.. Again, 
the value, of the phase difference «„, -«„ have been to a ^ [0;255J ^ Q 

representing and 255 representing 122^, aUovnng ^ values to be Splayed. The 

dark area, 7,0 and 71! have values close OT o. le representation, which 
corresponds to values for close to -„ ,„ m e phas e difference a.-^. whereas the 
light areas 720 and 72, have vah.es close to 255 i„ .he representation, which corresponds 

to values close to n in the phase difference n ry a j- , 

P ainerence or Bm - cc mk . Accordingly, a value of 128 in 

the representation corresponds to a value of 0 in the phase difference cc -„ The 
spread of values in the representation is shown in the histogram of Fig. 22B. The peaks 
730 and 73 1 in the histogram corresponds to phase difference a„ - a )Hk values. 

Step 150 determines a modulus of the phase difference [a^-a^y fig 23 A 
shows the modulus | a . m - ^ , of the phase dlfferencg represented in Fig. 22A 

Again the values have been adapted to a range [0;255] with 0 representing -n and 255 

127 

representing — n . As can be s«n from Fig. 23A, all the values have.a small variation. 
This is confirmed in the histogram shown in Fig. 23B, wherein the spread of values in the 
representation is shown. A sin gl e peak 910 in the histogram corresponds to modulus of 
Phase difference ,a m - a- , vaUles . „ is noted „„, ^ ^ ^ ^ ^ ^ ^ 

around the peak value. 

Step 660 calculates the estimated polarity factor H^, using Equahon (3?) 
The result shown in Fig. 24 wherem all regions have values -1 or 1, with the dark areas 
corresponding to values of -1 and the light areas corresponding to values of 1 . The phase 
differences y mn are calculated using Equation (59) and the relative phase-shift parameters 
"- |+ y -c-i> *» estimated in a consistent manner in step 670. 
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The relative phase-shift parameters S n are then substituted into the general PSA 
defined by Equation (39) and the fringe pattern offset «(,.,), ^ modulation b(x , y) , 
and phase are estimated in step 680. The resultant modulation b( x ,y), phase 

*W) and offset a( X>y ) of the fringe pattern/.^) of the example are shown in Figs 
5 25 A, 25B and 26 respectively. It is noted that Figs. 25A and 26 appear almost uniformly 
white, winch corresponds to a]most no variation in the values of the modulation *(*,„) 
and offset a( x ,y) with respect to coordinates for the example shown. 

Step 690 estimates the likely errors in the earlier contingent ana,ytic images g 
and tmk by means of a d Jre ct comparison of the analytic image ^ and the PSA image" 
• Step 692 determines whether these errors are higher than a preset threshold. If the errors 
^ smaller than the threshold, the method 600 ends in step 696. !f these errors 
than the threshold, the errors are incorporated into the weighting „ {x . y) m step o94 and 
the method 600 returns to step 670 in an iterative manner, where the new weighnng 

^ ^ aPPli6d t0 ^ CStimatl0n »~» Phase-difference Fmt , until the errors are 

sufficiently !ow. In a simplest embodiment, the weighting w[x>y) is set to , for _ 
ordznates ^yj where the comparison is good and 0 where the comparison is poor 

Fig. 27 shows an example where the comparison of step 690 yields a weighting 
WWChiS 1 ° Wh - th -~onisgoodand 0 .0where the comparison ls poor 

such p 15 generaUy n0t6d ** in EqUati ° nS ^ integralS ° V - a11 * and y 

such as Equates (41, 54, 55, 59 and 61), the integrals may be replaced by discrete 
summations for discretely sampled images. 

The implementation is applicable to interferon^ testing of optica! 
components and optical systems Th* „«,k a- "pucai 
t^u , , embodiment may also be used in non-destructive 

testing to study the deformation of objects under stress and the <t„H k • 
tn ;,w;a, u ■ ' the stud y of v »bratmg objects 

,de„„f y v lbranon mod es. Often, optical techniquM Me alM ^ visualise L 
flows, ^ variations in den5jty Mc jn ^ fomi ^ ^ patterns. 

ortenUiioT 7T " dMn ° dU,atine * P-em. of e Slimati „ 8 „ 

S ys.en, ,00, such as tha. shown i„ Fig . 28 wherejn > 
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may be implemented as software, such as an application program executing within the 
computer system 800. In particular, the steps of such methods are effected by instructions 
in the software that are carried out by the computer. The software may be stored in a 
computer readable medium, including the storage devices described below^- The software 
is loaded into the computer from the computer readable medium, and then executed by 
the computer. A computer readable medium having such software or computer program 
recorded on it is a computer program product. The use of the computer program product 
in the computer preferably effects an advantageous apparatus for the methods described 
above. 

The computer system 800 comprises a computer module 802, input devices such 
as a keyboard 810 and mouse 812, output devices including a printer 808 and a display 
device 804. 

The computer module 802 typically includes at least one processor unit 814, a 
memory unit 818, for example formed from semiconductor random access memory 
(RAM) and read only memory (ROM), input/output (I/O) interfaces including a video 
interface 822, and an I/O interface 816 for the keyboard 810 and mouse 812. A storage 
device 824 is provided and typically includes a hard disk drive 826 and a floppy disk 
drive 828. A magnetic tape drive (not illustrated) may also be used. A CD-ROM 
drive 820 is typically provided as a non-volatile source of data. The components 814 
to 828 of the computer module 802, typically communicate via an interconnected bus 830 
and in a manner which results in a conventional mode of operation of the computer 
system 100 known to those in the relevant art. Examples of computers on which the 
embodiments can be practised include IBM-PC's and compatibles, Sun Sparcstations or 
alike computer systems evolved therefrom. 

Typically, the application program of the preferred embodiment is resident on 
the hard disk drive 826 and read and controlled in its execution by the processor 814. 
Intermediate storage of the program may be accomplished using the semiconductor 
memory 818, possibly in concert with the hard disk drive 826. In some instances, the 
application program may be supplied to the user encoded on a CD-ROM or floppy disk 
and read via the corresponding drive 820 or 828, or alternatively may be read by the user 
from a network via a modem device (not illustrated). Still further, the software can also 
be loaded into the computer system 800 from other computer readable medium including 
magnetic tape, a ROM or integrated circuit, a magneto-optical disk, a radio or infra-red 
transmission channel between the computer module 802 and another device, a computer 
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readable card such as a PCMCIA card, and the Internet and Intranets including email 
transmissions and information recorded on websites and the like. The foregoing is merely 
exemplary of relevant computer readable mediums. Other computer readable mediums 
may be practiced without departing from the scope and spirit of the invention. 

The methods may alternatively be implemented in dedicated hardware such as 
one or more integrated circuits performing the functions or sub functions. Such dedicated 
hardware may include graphic processors, digital signal processors, or one or more 
microprocessors and associated memories. 

The foregoing describes only some embodiments of the present invention, and 
modifications and/or changes can be made thereto without departing from the scope and 
spirit of the invention, the embodiment being illustrative and not restrictive. 
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